Primary vs. Mets
Run Primary vs. Metastatic
#create the object
DESeq.herv.pm <- DESeqDataSetFromMatrix(countData = as.matrix(all_hervs[,rownames(pm_meta)]),
colData = pm_meta,
design = ~ Status)
## converting counts to integer mode
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeq.herv.pm$Status <- relevel(DESeq.herv.pm$Status, ref = "Primary")
Process count info
#size factor
DESeq.herv.pm <- estimateSizeFactors(DESeq.herv.pm)
#normalize & log-sized
DESeq.herv.pm_norm <- counts(DESeq.herv.pm, normalized = TRUE)
boxplot(log2(DESeq.herv.pm_norm+1), notch=FALSE, main = "Size-factor-normalized read counts", ylab="log2(read counts)", cex = .6, las = 2)
assay(DESeq.herv.pm,"log.norm.counts") <- log2(DESeq.herv.pm_norm+1)
#SD and PCA plot
msd.herv.pm <- vsn::meanSdPlot(DESeq.herv.pm@assays@data@listData$log.norm.counts, ranks=FALSE, plot = FALSE)
msd.herv.pm$gg + ggtitle("Sequencing depth normalized log2(read counts)") + ylab("standard deviation")
DESeq.herv.pm.trans <- DESeq(DESeq.herv.pm)
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 384 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
herv.pm.vst <- assay(varianceStabilizingTransformation(DESeq.herv.pm.trans))
herv.pm.pca <- pca(herv.pm.vst, metadata = colData(DESeq.herv.pm.trans), removeVar = 0.1)
## -- removing the lower 10% of variables based on variance
PCA Visualize
screeplot(herv.pm.pca, axisLabSize = 18, titleLabSize = 22) + ggtitle("SCREE Plot for HERVs, Primary + Metastatic")
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).
biplot(herv.pm.pca, lab = NULL, colby = "Cluster",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by Cluster")
biplot(herv.pm.pca, lab = NULL, colby = "Status",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by Status")
biplot(herv.pm.pca, lab = NULL, colby = "toMet",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by Metastasis Development")
biplot(herv.pm.pca, lab = NULL, colby = "Batch",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by Batch")
biplot(herv.pm.pca, lab = NULL, colby = "GNAQ_site", legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by GNAQ Mutation Site")
biplot(herv.pm.pca, lab = NULL, colby = "BAP1", legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by BAP1 Status")
biplot(herv.pm.pca, lab = NULL, colby = "BAP1_type",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by BAP1 Mutation Type")
biplot(herv.pm.pca, lab = NULL, colby = "Sex", legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by Sex")
biplot(herv.pm.pca, lab = NULL, colby = "SF3B1",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by SF3B1 Status")
options(ggrepel.max.overlaps = Inf)
biplot(herv.pm.pca, lab = herv.pm.pca$metadata$Sample, colby = "Cell_line",legendPosition = 'right') + ggtitle("PCA, HERVs Primary + Mets\nColor by Cell Line")
## Warning: ggrepel: 94 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
#ggsave("all HERVs Cell Line.png")
PCA Loading Plot
plotloadings(herv.pm.pca,
components = getComponents(herv.pm.pca, c(1,2)),
rangeRetain = 0.1,
labSize = 4.0,
title = 'Loadings plot for primary + metastatic HERVs',
subtitle = 'PC1, PC2',
caption = 'Top 1% variables',
shape = 24,
col = c('limegreen', 'black', 'red3'),
drawConnectors = TRUE)
## -- variables retained:
## MER101_16p12.2a, HERVH_13q22.3, HERVH_4p15.1a, PABLB_5q12.2, HERVH_10q23.32, HERVI_9q22.1, HERVH_6q24.1b, PRIMA41_13q22.3, MER4_Xq21.33c, HML6_3p22.3, ERVL_9q31.1, MER4B_9q31.1b, HERVL_21q22.2, ERVLE_2q24.3b, MER61_4q31.21b, HUERSP3B_Xq13.2b, ERV316A3_14q21.3c, HML5_12q12b, HERV9_6p21.2, HERVL_12q12b, HERVL_2q22.3a, HERVH_2q33.1a, MER41_15q13.3, HERVW_8q21.13, HERVIP10FH_5q32, MER4B_5q14.1, MER41_12q12, HUERSP3_6q22.31, MER41_5q15a, LTR25_11p14.1, ERVLB4_7q31.1c, ERVLB4_6q23.2a, ERVLE_15q12b, MER41_6q22.31b, HERVEA_4q22.1, MER61_16q12.1, HERVEA_Xp11.21, HERVFRD_5q11.2, HERVL_1q25.2, ERVLE_11q14.2c, HARLEQUIN_1q42.13a, MER41_15q21.2a, ERV316A3_14q21.3b, MER61_Xq11.1, HML3_7q36.3
#max.overlap = Inf
#DESeq.herv.pm$Cluster <- relevel(DESeq.herv.pm$Cluster, ref = "Metastatic")
DESeq.herv.pm <- DESeq(DESeq.herv.pm,parallel = T,test = "LRT", reduced = ~1)
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates: 6 workers
## mean-dispersion relationship
## final dispersion estimates, fitting model and testing: 6 workers
## -- replacing outliers and refitting for 401 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DE table on HERV
#resultsNames(DESeq.herv.pm)
#"Intercept" "Status_Primary_vs_Metastatic"
#Mets vs. All
herv_res_pm<-restable(min_fc = 2, max_p = 0.01, n_top_genes = 20, DESeq.herv.pm, cont = c("Status","Metastatic","Primary"))
## baseMean log2FoldChange lfcSE stat pvalue
## HML2_3q13.2 2.767446 5.240929 0.4065786 185.9204 2.470963e-42
## ERV316A3_19q13.2 12.567693 4.043621 0.3113588 179.5517 6.071500e-41
## HML2_11q23.3 8.365531 6.028128 0.4690936 178.1348 1.237872e-40
## HML2_19q11 20.850587 8.142469 0.5922696 174.7286 6.862310e-40
## HERVL_3q21.2 31.807864 5.196241 0.4022239 170.2177 6.631538e-39
## HERVIP10FH_17q21.32 6.563468 2.413798 0.1989844 162.2533 3.642169e-37
## HML2_3q27.2 8.151993 5.173398 0.4267370 159.7918 1.256427e-36
## HML3_7q36.3 183.126780 13.328096 0.7746528 148.7669 3.224660e-34
## HERVS71_11q13.2 29.115358 6.742577 0.5465766 142.4691 7.679163e-33
## ERVLE_1p34.3c 2.953931 4.629152 0.4146446 141.2407 1.425323e-32
## HML2_7q22.2 927.127296 3.950818 0.3492430 137.6423 8.725988e-32
## HML6_14q24.2 13.542162 5.113347 0.4463552 136.6914 1.408576e-31
## HERVI_16q22.2 79.982709 10.089581 0.7779156 134.5523 4.136857e-31
## ERVLE_17p11.2c 231.320369 -2.447170 0.1857741 133.0566 8.787103e-31
## LTR46_1p36.13 4.428574 5.966823 0.5536398 128.1465 1.042544e-29
## HERVL_21q22.2 10.671705 4.287555 0.3943519 128.1117 1.061013e-29
## PRIMA4_12p11.21a 7.244193 2.740667 0.2617433 121.6670 2.730203e-28
## MER4_11p15.4a 4.748947 2.391660 0.2390804 110.9362 6.110597e-26
## HML3_1p21.3 7.021049 3.324386 0.3386361 107.5918 3.302422e-25
## HERVE_6q15 5.955854 6.228608 0.6254671 107.1718 4.081993e-25
## padj
## HML2_3q13.2 9.300704e-39
## ERV316A3_19q13.2 1.142656e-37
## HML2_11q23.3 1.553116e-37
## HML2_19q11 6.457433e-37
## HERVL_3q21.2 4.992222e-36
## HERVIP10FH_17q21.32 2.284854e-34
## HML2_3q27.2 6.755989e-34
## HML3_7q36.3 1.517202e-31
## HERVS71_11q13.2 3.211597e-30
## ERVLE_1p34.3c 5.364916e-30
## HML2_7q22.2 2.985875e-29
## HML6_14q24.2 4.418234e-29
## HERVI_16q22.2 1.112224e-28
## ERVLE_17p11.2c 2.204977e-28
## LTR46_1p36.13 2.218697e-27
## HERVL_21q22.2 2.218697e-27
## PRIMA4_12p11.21a 5.408675e-26
## MER4_11p15.4a 1.150014e-23
## HML3_1p21.3 5.919198e-23
## HERVE_6q15 6.983918e-23
#Mets vs. All, excluding the batch hervs
herv_res_pm_filt <- restable(min_fc = 2, max_p = 0.01, n_top_genes = 20, DESeq.herv.pm, cont = c("Status","Metastatic","Primary"), exclude_gene = batch_hervs)
## Warning in if (exclude_gene != "") {: the condition has length > 1 and only the
## first element will be used
## baseMean log2FoldChange lfcSE stat pvalue
## HML2_3q13.2 2.767446 5.240929 0.4065786 185.9204 2.470963e-42
## ERV316A3_19q13.2 12.567693 4.043621 0.3113588 179.5517 6.071500e-41
## HML2_11q23.3 8.365531 6.028128 0.4690936 178.1348 1.237872e-40
## HML2_19q11 20.850587 8.142469 0.5922696 174.7286 6.862310e-40
## HERVL_3q21.2 31.807864 5.196241 0.4022239 170.2177 6.631538e-39
## HERVIP10FH_17q21.32 6.563468 2.413798 0.1989844 162.2533 3.642169e-37
## HML2_3q27.2 8.151993 5.173398 0.4267370 159.7918 1.256427e-36
## HML3_7q36.3 183.126780 13.328096 0.7746528 148.7669 3.224660e-34
## HERVS71_11q13.2 29.115358 6.742577 0.5465766 142.4691 7.679163e-33
## ERVLE_1p34.3c 2.953931 4.629152 0.4146446 141.2407 1.425323e-32
## HML2_7q22.2 927.127296 3.950818 0.3492430 137.6423 8.725988e-32
## HERVI_16q22.2 79.982709 10.089581 0.7779156 134.5523 4.136857e-31
## ERVLE_17p11.2c 231.320369 -2.447170 0.1857741 133.0566 8.787103e-31
## LTR46_1p36.13 4.428574 5.966823 0.5536398 128.1465 1.042544e-29
## HERVL_21q22.2 10.671705 4.287555 0.3943519 128.1117 1.061013e-29
## PRIMA4_12p11.21a 7.244193 2.740667 0.2617433 121.6670 2.730203e-28
## MER4_11p15.4a 4.748947 2.391660 0.2390804 110.9362 6.110597e-26
## HML3_1p21.3 7.021049 3.324386 0.3386361 107.5918 3.302422e-25
## HERVE_6q15 5.955854 6.228608 0.6254671 107.1718 4.081993e-25
## MER34B_1q23.3a 3.711454 4.097588 0.4150451 106.9418 4.584321e-25
## padj
## HML2_3q13.2 9.300704e-39
## ERV316A3_19q13.2 1.142656e-37
## HML2_11q23.3 1.553116e-37
## HML2_19q11 6.457433e-37
## HERVL_3q21.2 4.992222e-36
## HERVIP10FH_17q21.32 2.284854e-34
## HML2_3q27.2 6.755989e-34
## HML3_7q36.3 1.517202e-31
## HERVS71_11q13.2 3.211597e-30
## ERVLE_1p34.3c 5.364916e-30
## HML2_7q22.2 2.985875e-29
## HERVI_16q22.2 1.112224e-28
## ERVLE_17p11.2c 2.204977e-28
## LTR46_1p36.13 2.218697e-27
## HERVL_21q22.2 2.218697e-27
## PRIMA4_12p11.21a 5.408675e-26
## MER4_11p15.4a 1.150014e-23
## HML3_1p21.3 5.919198e-23
## HERVE_6q15 6.983918e-23
## MER34B_1q23.3a 7.502341e-23
print(paste0(nrow(herv_res_pm_filt), " of HERVs were reported to be significant."))
## [1] "845 of HERVs were reported to be significant."
write.csv(herv_res_pm_filt, "herv_res_pm_filt.csv")
QQ Plot
herv_res_pm_q <- results(DESeq.herv.pm) %>% as.data.frame()
Norm_quantile<- sort(-log10(seq(1, 1/nrow(herv_res_pm_q), length.out = nrow(herv_res_pm_q))), decreasing = TRUE)
herv_res_pm_q <- herv_res_pm_q[order(herv_res_pm_q$padj, decreasing = FALSE),]
herv_res_pm_q["Normal_Quantile"] <- Norm_quantile
ggplot(herv_res_pm_q, aes(x = Normal_Quantile, y =-log10(padj))) + geom_point()+geom_abline(intercept = 0, slope = 1, color="red") + xlab("Negative Log Normal Quantile (Expected)")+ylab("Negative Log 10 P Values (Observed)")+ggtitle("QQ Plot for HERV Primary vs Mets\nDesign=RECIST")
## Warning: Removed 2588 rows containing missing values (geom_point).
Heatmap
#pm_heat <- heattable(DESeq.herv.pm, herv_res_pm, select_string = "^(ERV|HERV|LTR|MER|HML).+",min_fc =2, max_p = 0.01)
pm_heat_filt <- heattable(DESeq.herv.pm, herv_res_pm_filt, select_string = "^(ERV|HERV|LTR|MER|HML).+",min_fc =2, max_p = 0.01,exclude_genes = batch_hervs)
pm_heat_clustering <- heattable(DESeq.herv.pm, herv_res_pm_filt, select_string = "^(ERV|HERV|LTR|MER|HML).+",min_fc =2, max_p = 0.01,select_genes = 50, exclude_genes = batch_hervs)
breaklist <- seq(-3, 3, by = 0.25)
#pheatmap(pm_heat, scale="row",main = "DE HERVs, All Primary vs. Pre-tx Metastatic\nCutoff p=0.01,Top 49 HERVs",fontsize = 20, size = 35, color = inferno(25), show_rownames = TRUE,cluster_cols = FALSE, cluster_rows = TRUE, breaks = breaklist, legend_labels = "Expression Z-score", annotation_col = pm_meta[,"Status", drop = FALSE], gaps_col = 80)
pheatmap(pm_heat_clustering, scale="row",main = "DE HERVs Clustering, Primary vs. Pre-tx Metastatic\nCutoff p= 0.01,Top 40 HERVs",fontsize = 20, size = 35, color = inferno(25), show_rownames = TRUE,cluster_cols = TRUE, cluster_rows = TRUE, breaks = breaklist, legend_labels = "Expression Z-score", annotation_col = pm_meta[,"Status", drop = FALSE], gaps_col = 80)
pheatmap(pm_heat_filt, scale="row",main = "DE HERVs, All Primary vs. Pre-tx Metastatic\nCutoff p= 0.01,Top 40 HERVs Excluding Batch",fontsize = 20, size = 35, color = inferno(25), show_rownames = TRUE,cluster_cols = FALSE, cluster_rows = TRUE, breaks = breaklist, legend_labels = "Expression Z-score", annotation_col = pm_meta[,"Status", drop = FALSE], gaps_col = 80)
png("DE HERVs Design Primary vs Metastatic.png", res = 1200, unit = "in", width = 12, height = 12)
pheatmap(pm_heat_filt, scale="row",main = "DE HERVs, All Primary vs. Pre-tx Metastatic\nCutoff p= 0.01,Top 40 HERVs Excluding Batch",fontsize = 20, size = 35, color = inferno(25), show_rownames = FALSE,cluster_cols = FALSE, cluster_rows = TRUE, breaks = breaklist, legend_labels = "Expression Z-score", annotation_col = pm_meta[,"Status", drop = FALSE], gaps_col = 80)
png("DE HERVs Design Primary vs Metastatic.Clustered.png", res = 1200, unit = "in", width = 12, height = 12)
pheatmap(pm_heat_clustering, scale="row",main = "DE HERVs Clustering, Primary vs. Pre-tx Metastatic\nCutoff p= 0.01,Top 40 HERVs",fontsize = 20, size = 35, color = inferno(25), show_rownames = FALSE,cluster_cols = TRUE, cluster_rows = TRUE, breaks = breaklist, legend_labels = "Expression Z-score", annotation_col = pm_meta[,"Status", drop = FALSE], gaps_col = 80)
Volcano
herv_res_all <- results(DESeq.herv.pm) %>% data.frame()
herv_res_all$hervs <- rownames(herv_res_all)
herv_res_all$hervs[(which(!(rownames(herv_res_all) %in% rownames(pm_heat_filt))))] <- ""
EnhancedVolcano(herv_res_all, lab = herv_res_all$hervs, x = 'log2FoldChange', y = 'padj', pCutoff = 0.01,FCcutoff = log2(2),title = "HERV Vocalno, Primary vs. Metastatic, p Cutoff = 0.01, Showing top 40 HERVs",titleLabSize = 30, labSize = 10,legendLabSize = 15,max.overlaps = Inf, drawConnectors = TRUE,arrowheads = FALSE,colAlpha = 1/2,col = c("black", "forestgreen", "yellow", "red2"),pointSize = 4)+theme( plot.title=element_text(size=45,face="bold"),axis.text=element_text(size=25),axis.title=element_text(size=25,face="bold"))
#ggsave("Valcano HERVs_PrimaryvsMets.png")
HERVs of Interest (17 Profile HERVs from Primary paper)
herv_profile <- c('ERVLE_17p11.2c','LTR46_Xq11.1','HERVH_6q24.1b','HERVH_3q21.3a','MER4B_3q21.2','HERVH_10q23.32','HML3_19q13.2','HERVI_9q22.1','MER4_1p21.3','HML3_3q21.2','ERVLE_5q13.2a','HERV3_19q13.42a','HERVE_Xp11.23','ERV316A3_6p21.33c','MER101_19q13.2c','ERVLE_2p13.2a','MER4B_12q22')
herv_res_17 <- herv_res_all[herv_profile,,drop = FALSE]
herv_res_17_sig <- subset(herv_res_17, padj <= 0.01)
print(paste0("There are ", nrow(herv_res_17_sig), " significantly DE HERVs between primary vs. Mets. As listed:" ))
## [1] "There are 9 significantly DE HERVs between primary vs. Mets. As listed:"
rownames(herv_res_17_sig)
## [1] "ERVLE_17p11.2c" "LTR46_Xq11.1" "HERVH_3q21.3a" "MER4B_3q21.2"
## [5] "HERVH_10q23.32" "HERVI_9q22.1" "ERVLE_5q13.2a" "ERVLE_2p13.2a"
## [9] "MER4B_12q22"
herv_res_17
## baseMean log2FoldChange lfcSE stat
## ERVLE_17p11.2c 231.320369 -2.44716980 0.1857741 133.05658977
## LTR46_Xq11.1 2084.308500 -1.59259318 0.1801353 65.81463790
## HERVH_6q24.1b 43.852405 -1.50449530 0.6871596 4.09095961
## HERVH_3q21.3a 2.978850 -1.99800006 0.5092303 13.76652668
## MER4B_3q21.2 74.287944 -0.99025872 0.2226502 17.96098383
## HERVH_10q23.32 127.738804 -3.60068196 0.4419191 44.21090281
## HML3_19q13.2 239.474843 -0.69242497 0.3171340 4.45184860
## HERVI_9q22.1 113.564864 -2.83310368 0.6061039 15.85432849
## MER4_1p21.3 56.698813 -0.63768465 0.2885318 4.61812195
## HML3_3q21.2 5.436287 -0.07509717 0.3636110 0.04305649
## ERVLE_5q13.2a 2.623870 1.42672476 0.2847621 27.90214455
## HERV3_19q13.42a 1.793846 0.73904242 0.3738196 4.61863662
## HERVE_Xp11.23 95.150772 0.27741746 0.2307454 1.49198307
## ERV316A3_6p21.33c 1553.189091 0.28266305 0.3919469 0.53288733
## MER101_19q13.2c 1.077808 1.21848724 0.7436845 3.33255678
## ERVLE_2p13.2a 3.594577 0.98863057 0.2876769 13.13696792
## MER4B_12q22 1.780568 2.51039744 0.4689431 33.06344866
## pvalue padj hervs
## ERVLE_17p11.2c 8.787103e-31 2.204977e-28 ERVLE_17p11.2c
## LTR46_Xq11.1 4.953894e-16 1.528398e-14
## HERVH_6q24.1b 4.311316e-02 6.991725e-02
## HERVH_3q21.3a 2.069917e-04 6.252943e-04
## MER4B_3q21.2 2.254795e-05 8.461664e-05
## HERVH_10q23.32 2.948346e-11 3.545551e-10
## HML3_19q13.2 3.486349e-02 5.814186e-02
## HERVI_9q22.1 6.840988e-05 2.335926e-04
## MER4_1p21.3 3.163586e-02 5.344586e-02
## HML3_3q21.2 8.356189e-01 9.101647e-01
## ERVLE_5q13.2a 1.276083e-07 7.435255e-07
## HERV3_19q13.42a 3.162637e-02 5.344586e-02
## HERVE_Xp11.23 2.219090e-01 2.941076e-01
## ERV316A3_6p21.33c 4.653955e-01 5.503119e-01
## MER101_19q13.2c 6.792121e-02 1.039672e-01
## ERVLE_2p13.2a 2.895260e-04 8.395808e-04
## MER4B_12q22 8.919992e-09 6.583304e-08
heat_17 <- heattable(DESeq.herv.pm, herv_res_17, select_string = ".+",min_fc = 0, max_p = 1,select_genes = 17, exclude_genes = "")
herv_res_17$Significance <- ifelse(herv_res_17$padj <= 0.01, "*","NS")
herv_res_17 <- herv_res_17[order(herv_res_17$padj),]
pheatmap(heat_17, scale="row",main = "Selected HERVs\nPrimary vs. Pre-tx Metastatic",fontsize = 20, size = 35, color = inferno(25), show_rownames = TRUE,cluster_cols = FALSE, cluster_rows = FALSE, breaks = breaklist, legend_labels = "Expression Z-score", annotation_col = pm_meta[,"Cluster", drop = FALSE], annotation_row =herv_res_17[,"Significance", drop = FALSE], gaps_col = 80)
Invidual HERV Plotting
Function to create error bars
data_sum <- function(df, variable_row, grouping_row){
g_factors <- as.factor(unique(df[,grouping_row]))
sum_table <- data.frame(matrix(nrow = length(g_factors), ncol = 2))
colnames(sum_table) <- c(variable_row, "sd")
rownames(sum_table) <- g_factors
for(i in 1:length(g_factors)){
sub_mean <- mean(df[which(df[,grouping_row] == g_factors[i]),variable_row])
sub_sd <- sd(df[which(df[,grouping_row] == g_factors[i]),variable_row])
sum_table[i,c(1:2)] <- c(sub_mean,sub_sd)
}
sum_table[grouping_row] <- g_factors
return(sum_table)
}
Significant count graphing
HML2_19q11 <- plotCounts(DESeq.herv.pm, gene = "HML2_19q11", intgroup = "Status", normalized = TRUE,returnData = TRUE)
HML2_summary<- data_sum(HML2_19q11, variable_row="count", grouping_row="Status")
ggplot(HML2_19q11, aes(x = Status, y = count)) +geom_point(aes(size = 1.5)) + ggtitle("Normalized Counts for HML2_19q11c\nSeparated by Status")+geom_errorbar(aes(x = Status, ymin = count-sd, ymax = count+sd),data = HML2_summary, width = 0.08, col = "grey50") + geom_errorbar(aes(x = Status, ymin = count, ymax = count),data = HML2_summary, width = 0.08, col = "red2") + theme( plot.title=element_text(size=16,face="bold"),axis.text=element_text(size=18),axis.title=element_text(size=20,face="bold"))
#ggsave("ERV316A3_3q13.12c Counts.png")
ERVLE_17p11.2c <- plotCounts(DESeq.herv.pm, gene = "ERVLE_17p11.2c", intgroup = "Status", normalized = TRUE,returnData = TRUE)
ERVLE_summary<- data_sum(ERVLE_17p11.2c, variable_row="count", grouping_row="Status")
ggplot(ERVLE_17p11.2c, aes(x = Status, y = count)) +geom_point(aes(size = 1.5)) + ggtitle("Normalized Counts for ERVLE_17p11.2c\nSeparated by Status")+geom_errorbar(aes(x = Status, ymin = count-sd, ymax = count+sd),data = ERVLE_summary, width = 0.08, col = "grey50") + geom_errorbar(aes(x = Status, ymin = count, ymax = count),data = ERVLE_summary, width = 0.08, col = "red2") + theme( plot.title=element_text(size=16,face="bold"),axis.text=element_text(size=18),axis.title=element_text(size=20,face="bold"))